library(Seurat)
library(tidyseurat)
library(harmony)
library(SingleR)
library(ggpubr)
library(ggrepel)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

key_cytokine <- c('Ccl20','Tslp','Flt3l','Csf2','Tnf','Il6')

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

# cavag24 -----
cavag24.path <- list.files('mission/FPP/psoria_imq/cavag24/', full.names = T)

cavag24.tb <- tibble(path = cavag24.path) |>
  mutate(type = str_extract(path, 'barc|feat|matr'),
         group = str_extract(path, '..(cont|IMQ|infe)')) |>
  pivot_wider(names_from = type, values_from = path) |>
  mutate(group = case_match(group,
                            '0_cont' ~ 'SA.ctrl',
                            '0_infe' ~ 'SA.infect',
                            '1_cont' ~ 'IMQ.ctrl',
                            .default = 'IMQ.psoria'))

cav24.mtx <- cavag24.tb |>
  pmap(function(group, barc, feat, matr){
    ReadMtx(mtx = matr, cells = barc,features = feat) |>
      add_name_field(group)}, .progress = T)

cav24.mtx[[1]][1:5,1:5]

glimpse(cav24.mtx)
cav24.merg <- RowMergeSparseMatrices(cav24.mtx[[1]], cav24.mtx[2:4])

sobj.cav24 <- cav24.merg |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj.cav24 <- sobj.cav24 |>
  PercentageFeatureSet('^mt-', col.name = 'mito.ratio')

sobj.cav24 |> VlnPlot('mito.ratio', pt.size = 0)

sobj.cav24 <- sobj.cav24 |>
  filter(mito.ratio < 5)

sobj.cav24 <- quick_process_seurat(sobj.cav24, leiden = F)

sobj.cav24 |> write_rds('mission/FPP/psoria_imq/cav24.rds')

mmus <- celldex::MouseRNAseqData()

sobj.cav24 <- sobj.cav24 |>
  mark_cell_type_singler(mmus, new_label = 'mmus_main')

sobj.cav24 |> DimPlot(group.by = 'mmus_main')

kc.marker <- c('Krt10','Krt14','Krt15','Krt16','Ivl','Dmkn','S100a14')

sobj.cav24 |> DotPlot(kc.marker, cluster.idents = T)

kc.fdot <- c(13,9,21,6,4,18,1,8,17,3)

sobj.cav24 |>
  mutate(kc.dotty = ifelse(seurat_clusters %in% kc.fdot,
                           'KC','nonKC')) |>
  DimPlot(group.by = 'kc.dotty')

kcm.de <- sobj.cav24 |>
  FindAllMarkers(features = kc.marker, only.pos = T, logfc.threshold = .5)

kcm.de |>
  as_tibble() |>
  filter(p_val_adj < .05) |>
  summarise(n(), .by = cluster) |>
  mutate(seurat_clusters = cluster, kc.5score = `n()`, .keep = 'none') |>
  left_join(x = sobj.cav24, y = _) |>
  FeaturePlot('kc.5score',cols = c('blue','red'))

sobj.cav24 <- sobj.cav24 |>
  mutate(cell.type = ifelse(seurat_clusters %in% kc.fdot,
                           'Keratinocytes', mmus_main))

sobj.cav24 |>
  DimPlot(group.by = 'cell.type', cols = 'Paired')

sobj.c24kc <- sobj.cav24 |>
  filter(cell.type == 'Keratinocytes')

sobj.c24kc <- sobj.c24kc |>
  quick_process_seurat(leiden = F, skip_norm = T)

sobj.c24kc |>
  DotPlot(kc.marker, cluster.idents = TRUE)

sobj.c24kc |>
  DotPlot('Trpv3', cluster.idents = TRUE)

FindAllMarkers(sobj.c24kc, features = 'Trpv3', only.pos = T) |>
  filter(p_val < .05)

v3.fc.ctrl <- sobj.c24kc |>
  filter(str_detect(orig.ident, 'ctrl')) |>
  FindAllMarkers(features = 'Trpv3', logfc.threshold = 0)

v3.fc.ctrl |>
  ggplot(aes(cluster, avg_log2FC)) + geom_col()

FeaturePlot(sobj.c24kc, features = 'Trpv3')

sobj.c24kc <- sobj.c24kc |>
  mutate(cell.type = ifelse(seurat_clusters %in% c(1,2,4,13),
                            'Trpv3-hi-keratinocytes', 'Trpv3-lo-keratinocytes'))

sobj.c24kc |> DimPlot(group.by = 'cell.type')

sobj.cav24 |> filter(str_detect(orig.ident, 'IMQ')) |>
  DimPlot(group.by = 'orig.ident')

sobj.c24kc
## Key cytokine ------
sobj.c24kc |> VlnPlot('Banf1', group.by = 'orig.ident', pt.size = 0)

sobj.c24kc |> DotPlot(c('Banf1', key_cytokine))

Idents(sobj.c24kc) <- 'cell.type'

kcv3.list <- sobj.c24kc$cell.type |> unique()

kcv3.list <- kcv3.list |>
  set_names(kcv3.list)

imq.deg.c24 <- kcv3.list |>
  map(\(x)FindMarkers(sobj.c24kc, subset.ident = x, group.by = 'orig.ident',
              features = c(key_cytokine, 'Gata6'),
              ident.1 = 'IMQ.psoria', ident.2 = 'SA.ctrl') |>
        as_tibble(rownames = 'gene')) |>
  list_rbind(names_to = 'cell.type')

imq.deg.c24 |>
  ggplot(aes(gene, avg_log2FC, fill = p_val_adj < .05)) +
  geom_col() +
  facet_wrap(~cell.type) +
  theme_pubr() +
  scale_fill_manual(values = c('grey','red')) +
  labs(title = 'Key cytokine differential expression in KC after IMQ psoriasis',
       fill = 'Significant?',
       subtitle = 'GSE230512')
  
### TF regulating Ccl20? ------
sobj.c24kc |> VlnPlot('Banf1', group.by = 'cell.type')

sobj.c24kc |> FindAllMarkers(features = 'Banf1')

## MVA --------
sobj.c24kc |> DotPlot(kegg_mva)

imq.deg.c24 <- kcv3.list |>
  map(\(x)FindMarkers(sobj.c24kc, subset.ident = x, group.by = 'orig.ident',
                      features = kegg_mva[kegg_mva != 'Idi2'],
                      ident.1 = 'IMQ.psoria', ident.2 = 'SA.ctrl') |>
        as_tibble(rownames = 'gene')) |>
  list_rbind(names_to = 'cell.type')

imq.deg.c24 |>
  left_join(mva_fct) |>
  ggplot(aes(ordered, cell.type, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  theme_pubr(legend = 'right', x.text.angle = 45) +
  scale_color_distiller(palette = 'RdBu', limits = c(-7.5,7.5)) +
  labs(title = 'Key cytokine differential expression in KC after IMQ psoriasis',
       fill = 'Significant?',
       subtitle = 'GSE230512')

# qu23 -----
path.q23 <-
list.files('mission/FPP/psoria_imq/', full.names = T, recursive = T) |>
  str_subset('qu23')

tb.q23 <- tibble(path = path.q23) |>
  mutate(type = str_extract(path, 'barc|gene|matr'),
         group = str_extract(path, '(SBM|IMQ)')) |>
  pivot_wider(names_from = type, values_from = path)
  
mtx.q23 <- tb.q23 |>
  pmap(function(group, barc, gene, matr){
    ReadMtx(mtx = matr, cells = barc,features = gene) |>
      add_name_field(group)}, .progress = T)

mtx.q23[[1]][1:5,1:5]

glimpse(mtx.q23)
merg.q23 <- RowMergeSparseMatrices(mtx.q23[[1]], mtx.q23[[2]])

sobj.q23 <- merg.q23 |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

library(magrittr)
sobj.q23 %<>% PercentageFeatureSet('^mt-', col.name = 'mito.ratio')

sobj.q23 |> VlnPlot('mito.ratio', pt.size = 0)

sobj.q23 %<>% filter(mito.ratio < 5)

sobj.q23 %<>% quick_process_seurat(leiden = F)

sobj.q23 |> DimPlot(group.by = 'orig.ident')

sobj.q23 %<>% mark_cell_type_singler(mmus, new_label = 'mmus_main')

mmus$label.fine |> table()

sobj.q23 |> DimPlot(group.by = 'mmus_main')

sobj.q23 |> DotPlot(kc.marker, cluster.idents = T)

sobj.q23 %<>%
  mutate(cell.type = ifelse(seurat_clusters %in% c(3,4,13,15),
                            'Keratinocytes', mmus_main))

sobj.q23 |> DimPlot(group.by = 'cell.type')

sobj.q23 |> write_rds('mission/FPP/psoria_imq/q23.rds')

kc.q23 <- sobj.q23 |> filter(str_detect(cell.type, 'Kera'))

kc.q23 %<>% quick_process_seurat(leiden = F, skip_norm = T)

kc.q23 |> DotPlot('Trpv3', cluster.idents = T)

kc.q23 |> FindAllMarkers(features = 'Trpv3', only.pos = T)

kc.q23 %<>% mutate(cell.type = ifelse(seurat_clusters %in% c(1,7),
                                      'Trpv3-hi-KC','Trpv3-lo-KC'))

kc.q23 |> DimPlot(group.by = 'cell.type')

## cytokine -----
kc.q23 %<>% mutate(subgroup = str_c(cell.type, '_', orig.ident),
                   cell.type = fct_relevel(cell.type, 'Trpv3-lo-KC'))

kc.q23 |> DotPlot(key_cytokine, group.by = 'subgroup')

kc.q23 |> filter(orig.ident == 'IMQ') |>
  VlnPlot(c('Ccl20','Il6'), group.by = 'cell.type')

kc.q23 |>
  FindMarkers(features = key_cytokine[key_cytokine != 'Flt3l'],
              group.by = 'cell.type',
              ident.1 = 'Trpv3-hi-KC')

## MVA ------
kc.q23 |> DotPlot(kegg_mva, group.by = 'subgroup') +
  rotate_x_text(45) +
  ggtitle('MVA expression in KC')

## TF -------
kc.q23 |> DotPlot(c('Gata6','Banf1','Rcor1'), group.by = 'subgroup')

# qu24 ---------
path.q24 <-
  list.files('mission/FPP/psoria_imq/', full.names = T, recursive = T) |>
  str_subset('qu24')

tb.q24 <- tibble(path = path.q24) |>
  mutate(type = str_extract(path, 'barc|gene|matr'),
         group = str_extract(path, '(CE|IMQ|NORMAL)')) |>
  pivot_wider(names_from = type, values_from = path)

mtx.q24 <- tb.q24 |>
  pmap(function(group, barc, gene, matr){
    ReadMtx(mtx = matr, cells = barc,features = gene) |>
      add_name_field(group)}, .progress = T)

mtx.q24[[1]][1:5,1:5]

glimpse(mtx.q24)
merg.q24 <- RowMergeSparseMatrices(mtx.q24[[1]], mtx.q24[-1])

sobj.q24 <- merg.q24 |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj.q24 %<>% PercentageFeatureSet('^mt-', col.name = 'mito.ratio')

sobj.q24 |> VlnPlot('mito.ratio', pt.size = 0)

sobj.q24 %<>% filter(mito.ratio < 5)

sobj.q24 %<>% quick_process_seurat(leiden = F)

sobj.q24 %<>% mark_cell_type_singler(mmus, new_label = 'mmus_main')

sobj.q24 |> DotPlot(kc.marker, cluster.idents = T)

kc.dotty = c(19,13,30,15,4,3,5,6,18)

sobj.q24 %<>% mutate(cell.type = ifelse(seurat_clusters %in% kc.dotty,
                                      'Keratinocytes', mmus_main))

sobj.q24 |> DimPlot(group.by = 'cell.type')

sobj.q24 |> write_rds('mission/FPP/psoria_imq/q24.rds')

sobj.q24 <- read_rds('mission/FPP/psoria_imq/q24.rds')

## Kera ------
kc.q24 <- sobj.q24 |>
  filter(str_starts(cell.type, 'Kera'))

kc.q24 %<>% quick_process_seurat()

kc.q24 |> DotPlot(kc.marker, cluster.idents = T)

kc.q24 |> DotPlot('Trpv3', cluster.idents = T)

kc.q24 |> FindAllMarkers(features = 'Trpv3', only.pos = T)

kc.q24 %<>% mutate(cell.type = ifelse(seurat_clusters %in% c(3),
                                      'Trpv3-hi-KC','Trpv3-lo-KC'))

kc.q24 |> DimPlot(group.by = 'cell.type')

kc.q24 %<>% mutate(subgroup = str_c(cell.type, '_', orig.ident),
                   cell.type = fct_relevel(cell.type, 'Trpv3-lo-KC'))

kc.q24 |> write_rds('mission/FPP/psoria_imq/q24.kc.rds')

## cytokine -------
kc.q24 |> 
  filter(!str_detect(subgroup, 'CE')) |>
  DotPlot(key_cytokine, group.by = 'subgroup', cols = 'RdYlBu')

## MVA -----
kc.q24 |>
  filter(!str_detect(subgroup, 'CE')) |>
  mutate(subgroup = str_replace(subgroup, 'NORMAL', 'ctrl')) |>
  DotPlot(kegg_mva, group.by = 'subgroup', cols = 'RdYlBu', dot.scale = 2) +
  labs(title = 'MVA expression in KC', x = 'Gene', y = 'Cell type') +
  theme_jpub +
  rotate_x_text(45)

publish_pdf('mission/FPP/psoria_imq/imq.kc.mva.bubble.pdf',
            width = 65, height = 40)

## TF ---
kc.q23 |> DotPlot(c('Gata6','Banf1','Rcor1'), group.by = 'subgroup')

library(data.table)

cr.frag71 <- fread('~/append-ssd/alaria2/sun2024psoria/pso_ctrl_out/outs/fragments.tsv.gz')

str1 <- c("ABCDEF", "GHIJKL", "MNOPQR")
str2 <- c("ABCGEF", "GHIJKM", "MNOPQS", "ABCDEF", "GHIJKL", "MNOPQR")

# Find strings in str2 with max Hamming distance of 1 from str1
result <- str2[apply(adist(str1, str2, costs = c(ins = 0, del = 0, sub = 1)), 1, function(x) any(x == 1))]

print(result)

str2[unlist(lapply(str1, function(x) which(adist(x, str2, costs = c(ins = 0, del = 0, sub = 1)) == 1)))]

adist(str1, str2)
hamming6 |> rowMin()
hamming1k <- adist(head(cmp.frag71$V4, 1000), cratac$bc) |> rowMin()
